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Abstract. In recent years, quantities derived from the heat equation 
have become popular in shape processing and analysis of triangulated 
surfaces. Such measures are often robust with respect to different kinds 
of perturbations, including near-isometries, topological noise and partial- 
ities. Here, wc propose to exploit the semigroup of a Schrodinger operator 
in order to deal with texture data, while maintaining the desirable prop- 
erties of the heat kernel. We define a family of Schrodinger diffusion 
distances analogous to the ones associated to the heat kernels, and show 
that they are continuous under perturbations of the data. As an appli- 
cation, we introduce a method for retrieval of textured shapes through 
comparison of Schrodinger diffusion distance histograms with the earth's 
mover distance, and present some numerical experiments showing supe- 
rior performance compared to an analogous method that ignores the 
texture. 
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1 Introduction 

There is an ever growing quantity of 3D shapes available, either scanned from 
real objects, manually modelled by artists, or acquired from other sources. Ad- 
equately classifying them, and being able to find similar and dissimilar models 
is therefore increasingly important, and automatic solutions are needed for the 
goal of efficient computerized shape retrieval. 

For a retrieval method to be useful, given the variability of shapes, often 
some invariancc properties are required. The most obvious one is translation 
and rotation (that is, Euclidean) invariancc. Their scale is often arbitrary, so 
sometimes it is also interesting to enforce invariance with respect to global or 
local scaling. More challenging is recognition in classes of non-rigid shapes, like 
shapes representing human faces, animals or animated characters. In these cases, 
only the intrinsic geometry can be used, thus enforcing invariance to isometries, 
and robustness with respect to near-isometries. 

Often, geometric models include textures, which arc an intc!gral part of the 
representation. It is therefore natural to try and use the texture information to 
better distinguish between objects, for example in cases like separating between 
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a horse and a zebra, classifying different species of fish with similar shapes but 
different colors and patterns, or categorizing archaeological findings. 

In this paper, we introduce a representation that incorporates texture data 
within several recent methods of shape analysis and retrieval, which them- 
selves depend only on intrinsic geometry and hence are appropriate for non- 
rigid shapes. Our method fully inherits the desirable invariance and robustness 
properties of these methods, while also utilizing the texture of the shapes. 

The paper is organized as follows. In Section 2 we briefly review previous 
efforts and basic concepts on which our method is based. In Section 3, wc define 
our central quantities, a family of diffusion distances based on diffusion with 
Schrodinger operators incorporating the texture data, and present some theo- 
retical results about them. Then, in Section 4, wc present a system of shape 
retrieval based on comparison of histograms of Schrodinger diffusion distances 
with the earth mover's distance. Finally, in Section 5, we present some experi- 
mental results obtained with our representation model. 

2 Diffusion in Shape Analysis. Previous Works 

Adding to a long history of use of Laplace operators in geometry processing 
applications [1], the spectral decomposition of Laplace-Beltrami operators on 
surfaces has proven to be useful for tasks of shape analysis and comparison [2]. 

In [3] , diffusion through the heat equation, constructed from the spectral de- 
composition, was used for comparison of shapes through the introduction of the 
heat kernel signature (HKS). Since then, many methods have used descriptors 
for shapes built from the heat kernel [4] [5] . 

Diffusion distances were introduced by Coifman and Lafon in [6] for data 
analysis of point clouds, under the basic assumption that the sampled points 
come from an underlying low-dimensional manifold. Recently, diffusion distances 
have also received considerable attention for shape analysis and retrieval tasks, 
for example in shape recognition [7] or shape matching [8]. 

Recently, an approach to shape retrieval including texture data was proposed 
in [9], introducing three channels of texture (in the Lab color space) through 
a higher dimensional embedding, similar to the Beltrami framework [10]. In 
comparison, the method presented here supports a single channel for the texture, 
yet it has a clear interpretation in terms of diffusion on the original shape. It also 
requires lower order derivatives of the texture (at most one in our case, versus 
two for the embedding approach), making it less sensitive to noise. 

3 Schrodinger Operators and Diffusion 

We will consider our surfaces to be compact two-dimensional manifolds. denotc;d 
by M and embedded in R^, with triangular meshes as discretizations. In terms 
of a local parametrization and the corresponding first fundamental form g, one 
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can define the Laplace-Beltrami operator through the formula 

A,f = (div, o grad,)(/) = J2 -jr=A {V\9\9''djf) , (1) 



^,J V 151 

where \g\ is the determinant of the metric, g^^ are the components of the inverse 
of the metric, and di denotes partial derivative with respect to the i-th coordi- 
nate. It is well known [11] that the Laplace-Beltrami opcirator doesn't depend 
on the coordinate functions chosen, and since it's defined in terms of g, it is in- 
variant under transformations that preserve g, that is, isometries. This operator 
is a generalization of the standard Laplacian in M", for many of the processes 
associated to the former, like difi^usion and smoothing. 

Consider a function F : M — )• K, which we require to be bounded, but without 
needing any further regularity, in particular not necessarily differentiable or even 
continuous, to which we will refer as the potential. A Schodinger operator on the 
surface M is an operator of the form Ag — V, with V being considered as a 
multiplication operator, that is, Vf{x) = V{x)f{x). 

These operators are referred to as Schrodinger operators, and play a major 
role in quantum mechanics, where the study of their spectrum is key to un- 
derstanding the Schrodinger equation. In those cases the potential is usually 
unbounded, which makes the analysis a challenge of its own. One can also con- 
sider the diffusion equation associated to these operators, 

idtu{x, t) = Agu{x, t) - V{x)u{x, t) 
|u(a:;,0) = Uq{x), 

which will be the equation defining the quantities that we will use in what 
follows. Boundary conditions are not needed since M is compact. The following 

result asserts that for these operators, on the continuous level, everything works 
as expected, mimicking the situation with the Laplacian and heat kernels: 

Theorem 1 Let {M,g) be a compact Riemannian manifold of class C^, Ag the 
Laplace-Beltrami operator on (M, g), and V G L°°{M, fig), where fXg is the mea- 
sure associated to the Riemannian volume element, with V >Q. Then, the oper- 
ator —Ag + V admits a spectral decomposition {{(l>i,^j)}'^i, such that {4>j}'jLi 
is an orthonormal basis for L'^{M, jig), Xj > and lim = +00. 

j-i-oo 

Moreover, there exists a family of functions ht G L'^{M,ijg ), such that for all 
uo e Z/^(M, yUg), the unique solution of (2) is given by 

u{x,t) = / ht{x,y)uo{y)diJ,g{y), (3) 
Js 

and the following formula holds 

00 

ht{x,y)=Y,e-^''Mx)My)- (4) 

i=l 



4 



Jose A. Iglesias, Ron Kimmel 



Proof. The proof is essentially the same as in the Laplacian case. We refer to [12] 
for a standard proof. Let us note, however, that both the hypotheses that the 
manifold is compact and the potential V is bounded are essential, as otherwise, 
the discreteness of the spectrum is not guaranteed, since in those cases the 
involved resolvents could fail to be compact. 

Based on this result, a squared diffusion distance [6] can be defined as the 
norm of the difference of the kernels for Equation (2), and using Equation (4) 

oo 

4{x,y) = \\ht{x,-) - h{y,.)\\l. =Y,e-^^^\U^) ' Uv)?- (5) 

j=i 

Note, that the sign in the exponential arises because wc have defined Aj > 
to be the eigenvalues of —{Ag — V). These distances enjoy the same properties 
as those associated to the heat kernel, since all the properties proved in [6], 
including the fact that the formula above defines a distance, are valid for more 
general semigroups and not just the one associated to the Laplacian. 

We also have that the solutions to the diffusion equation, and therefore any 
quantities derived from it, in particular our diffusion distances, are continuous 
with respect to perturbations of the potential. Namely, 

Proposition 2 Consider the problem: 

dtu{x, t) ^ Au{x, t) - {V{x) + eN{x))u{x, t) 

.(x,0)=.o(x), 

where V, N e L°°{M,^ig), V >Q and V + eN > {] for some e > 0. Then, the 
solutions to of the above problem converge strongly in L^{M,Hg) to the ones of 
problem 2 as e ^ 0, for each fixed t > 0. 

Proposition 2 can be proved by using standard results in perturbation theory 
of linear semigroups. We provide a proof in appendix A. 



3.1 The Feynman-Kac Formula 

To shed some light on the behavior of the solutions to Equation (2), we provide an 
informal discussion on a well-known stochastic interpretation of such solutions. 

The Feynman-Kac formula [13] expresses the solution of a diffusion equation 
in terms of Brownian motion (strictly, the Wiener process on our space, X) 

uix, t)^E (^wo(^o) exp ^ V{Xr )dT^ \Xt=x^, (7) 

the conditional expectation meaning that we take averages of all the Wiener 
paths that reach x at time t, starting from elsewhere. Note that the integral 
inside the exponential involves the Wiener process itself, and hence needs to be 
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understood in the sense of stochastic integrals. A rigorous treatment of this is 
beyond the scope of this paper, but can be found in [13]. 

In the case of the heat equation, ^ = 0, the initial values are transported 
over random paths, and the expected value over all paths that reach a point at 
a given time is the value of our solution. This kind of averaging property is the 
reason behind the robustness to different kinds of noise that diffusion distances 
and other quantities derived from heat kernels possess. 

For Schrodinger operators, one can think of this transported value being 
modulated exponentially by the potential V{x) > 0, in a way consistent with 
what one gets by disregarding the diffusion term to end up with ut + Vu — 0. 
The transported values will be decreased according to how large V is in the areas 
that the Brownian motion crosses, on average. Figure 1 illustrates this behavior. 




Fig. 1. From left to right, row-wise: Potential and kernel ht{p,-) at three different 
times. 



4 Textured Shape Retrieval with Schrodinger Diffusion 

In our context, all references to textures on shapes will in fact be about vertex 
colorings. This is because only at that level the notions of a mapping on the 
surface and the diffusion can have discretizations consistent with one another, 
as required for our Schrodinger operators. In the case of textures mapped on 
triangles, one could induce a vertex map by averaging over the Voronoi region 
corresponding to the vertex in its one-ring neighborhood, for example. 

For our shape retrieval application, we intend to define a distance between 
signatures of the shapes, which in our case will be histograms of Schrodinger 
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diffusion distances between points. The retrieval method woiild then select from 
a database the shapes with the smallest distance to the query. In what follows, 
we assume our texture is a diffcrcntiable function / : M M in the continuous 
model, and a vertex function I{vi) for the discretized version. 

4.1 Operator Discretization 

A popular discretization for the Laplace-Beltrami operator of a surface is the so- 
called cotangent weight scheme [14] , where the weights arc the sums of cotangents 
of angles adjacent to the edge, and the normalization coefhcients are the Voronoi 
areas corresponding to the vertex. One of its disadvantages is the fact that 
for meshes with obtuse angles, the edge weights become negative, making it 
unsuitable [15] for simulating diffusion processes. 

So instead, we chose the 'Mesh Laplacian' discretization of [16], inspired by 
the one introduced by Belkin and Niyogi for data analysis of point clouds in [17]. 
For triangular meshes, it is given by 

TeT{Vi) wgt 

where T(wi) denotes the set of triangles in a neighborhood of Vi, A{t) is the area 
of the triangle r and s is a scaling parameter. 

This discretization has several advantages. The ones that are the most useful 
for us are that the weights are nonnegativc by definition, and the fact that it con- 
verges pointwise to the continuous Laplace-Beltrami operator, when the meshes 
approximate a smooth surface and T{vi) is always the whole shape, as proved 
in [16]. Spectral convergence in a probabilistic sense, of the point-cloud version 
of the operator, is proved in [18]. In our case, since it would not be practical to 
use the whole shape as a neighborhood, we have taken the neighborhoods Tivi) 
to be the one-ring neighborhood of each vertex. The parameter s was choosen in 
a uniform way, not taking into account the size of the different neighborhoods, 
as one fifth of the median of the edge lengths over the whole shape. 

After discretizing the Laplace-Beltrami operator as above, our discrete oper- 
ators are defined in the obvious way, 

{{A, - V)f){vi) = AJivi) - V{vi)f{vi). (9) 

4.2 Choice of Potential and its Discretization 

One choice for the potential V would be just to take V = I, I corresponding 
for example to the luminance of the texture. This would have the advantage of 
not having to explicitly compute any derivative of /, but in turn would make it 
depend on the reference taken for the texture, that is, on transformations of the 
kind I = I + c, where c is a constant. 

Another option is to use an edge indicator for the textures as the potential V, 
the most straightforward being the modulus of the surface gradient, V = | V/|, V 
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being the (Ricmannian) gradient on the surface. Other options are V = log(l+/) 
or y = log(l + |V/|), as a way to mitigate the exponential decay caused by the 
potential V in Formula (7). In our experiments below, we used this last potential. 

Intuitively, from the random walk interpretation and the Feynman-Kac For- 
mula (7), we see that it will be harder to diffuse across edges of the texture, 
while in constant areas, the behavior will be that of the usual heat equation. 
This can be appreciated in Figure 2. 

We will now describe the approximation of the gradient of a function on the 
surface that we employed, which is the same as in [19]. Let v be the vertex we 
are interested in, and r = {u, v, w} a triangle of the shape having ii as a vertex. 
Denote 



U — V W — V 

\u — v\ \w — v\ 



T 

r 

T 

r 



(10) 



where r^, denote transposes. Note, that P is the change of basis matrix from 
the canonical basis to ru,rw, which form a basis for the plane that contains 
the triangle. Also, denote by the values of our function on the three 

vertices. Then, we can consistently approximate the norm of the gradient by 

(iv/i:)^ = f^) {pp^r (m) ■ (11) 

Finally, to obtain our discretization, we average this approximation over the 
one-ring neighborhood of v, 



4.3 Point Selection and histogram comparison 

To obtain a signature based on diffusion distances, one needs to select points 
between which the distances are computed. In our case, we performed farthest 
point sampling [20] based on Euclidean ambient distances (for simplicity) to pick 
100 points, for which we computed the diffusion distance map to the rest of the 
shape. Then, the results for each of them are combined in a global histogram, 
which will serve as descriptor for the shape. The histograms were quantized with 
120 bins, and normalized, to compare shapes with any number of vertices. 

After obtaining the histograms of diffusion distances to be used as signatures, 
we need a way to compare them consistently. A popular approach for comparing 
probability distributions (and hence normalized histograms) is the Earth Mover's 
Distance, or EMD [21]. We use it to compute distances between the histograms of 
diffusion distances, as a means of comparing them. Let us note that even if in the 
one-dimensional continuous case computing the optimal cost is straightforward 
from the cumulative distribution functions, in the discrete case one still needs to 
solve a flow network optimization problem, which is computationally expensive. 
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Fig. 2. Diffusion on a textured shape, witli V — a{l + /?j V/j). Left to right, row wise: 
textured shape, modulus of the gradient, snapshot of kernel, diffusion distance. Observe 
that some of the stripes of the zebra can be seen in the distance and the kernel is shaped 
by them, since texture edges work against the diffusion. The source is marked in white. 



if done naively [22]. In our case, we have used the method and code of Pele and 
Werman [23] to efficiently approximate the EMD between our histograms. 

A similar approach to shape retrieval, but with local descriptors, was used in 
[24], where the distances used were the inner distances inside planar shapes. A 
more sophisticated approach using histograms of geodesic distances and Wasser- 
stein metrics can be found in [25]. We emphasize that we purposefully used a 
classification method which is far from being state of the art, to be able to better 
demonstrate the increased performance with just a moderate-sized database. 



5 Experimental Results 

For our experiments, the shapes were taken from the TOSCA nonrigid shape 
database [26], and the textures were manually added as vertex colorings. The 
final database consisted of 73 shapes, belonging to 8 different classes. Inside each 
class, the shapes differed by an almost-isometric deformation (different 'poses'). 

The Laplace-Beltrami operators were discretized through the scheme de- 
scribed above, before combining them with the norms of the gradient of the 
texture data, as described in Section 4. Diffusion distances were computed di- 
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rectly from the definition, using 100 eigenvalues and eigenf unctions, which were 
in turn computed from the corresponding matrices of the discrete operators. 

Then, the EMD between the histograms from each pair of shapes was com- 
puted, and finally the two closest matches for each shape were chosen as candi- 
dates for retrieval. This whole process was done for V" = and V — alog(l + 
^|V/|), with a and /3 normalization constants which control the resistance to 
diffusion induced by the texture. 

Results for the distances are shown in Figure 3, and Table 1 shows the amount 
of correct matches and averages of distances inside and outside the classes, after 
normalization with respect to the maximum distance between elements. 

As we can see, the average distance between shapes in the same class is 
reduced, while the average distance between shapes in different classes increases, 
providing better separation between the classes. Also, some queries that would 
produce incorrect matches with just the geometry are correctly matched when 
also using the texture. 




Fig. 3. Distance matrices between the signatures of the shapes. Lines indicate separa- 
tion between classes. White corresponds to zero distance, black to maximum distance. 
Left: Results without texture. Right: Results with texture. Lower: Representatives of 
each class of shapes, in the same order as the matrices. 



Table 1. Numerical results for both cases 





= 


T/ = alog(l-|-/3|V7|) 


Nearest shape belongs to correct class 


52/73 


65/73 


Second nearest shape belongs to correct class 


41/73 


50/73 


Normalized avg. distance for shapes in same class 


0.1440 


0.1167 


Normalized avg. distance for shapes in different classes 


0.2905 


0.3822 
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6 Conclusions 

We defined a family of diffusion distances based on the diffusion associated to 
Scfirodinger operators, for use in triangulated shapes with textures, where quan- 
tities derived from the texture are introduced as the potential part of the oper- 
ator. These are at least continuous with respect to perturbations of the texture. 

The practical usefulness of our method was illustrated by a simple retrieval 
example using global diffusion distance histograms as descriptors. Using the 
available texture information resulted in better performance than using only the 
geometry data through standard heat diffusion. 

Our approach could also be useful to incorporate texture data in other meth- 
ods of shape analysis using Laplacian operators but not heat diffusion explicitly, 
such as the ones in [27] and [2]. 

Acknowledgments. This work was supported by the European Commission 
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A Proof of Proposition 2 

We want to apply [28] IX. 2. 16, which tells us that for strong convergence of a 
semigroup e*-^" — >■ e*-^ as e — 0, it is enough that for some ^ > larger than the 
exponential growth of the semigroup, the corresponding resolvents (T^ — 
converge strongly. In this case, since our operators T = Ag — V {V ^ {M, jig)) 
are m-accretive (the semigroup is contractive), any ^ > will suffice. Choose 
/ e L^{M,Hg) with ||/||i2 = 1, C > + ||iV||Loo such that £, is not an 

eigenvalue of —Ag + V, and consider the equations 



Note, that by the Fredholm alternative ([29], 8.95) and the fact that the spectrum 
of —Ag + {V + eN) is nonnegative, there is a unique solution for both equations. 
We can also assume e < 1. Multiplying the first equation of (13) by and 
integrating, we obtain 



460-467 



-AgU, + {V + eN)u, + ^u, = f 
—AgW + Vw + = f. 



(13) 




12 



Jose A. Iglesias, Ron Kimmel 



so that 




ll/IUHkelU^ 



< 



Me II 1,2 



(14) 



C-e||7V||^oc-||F|| 



ll^ll 



and hence is bounded in L'^{M,fig) independently of e. Without loss of gen- 
erality we can consider e = 1/n, and by elliptic regularity, we can take a subse- 
quence Uk such that ?i„ j. converges weakly to some uq G H'^{M). Now substract 
both equations in (13), multiply by an arbitrary h in H^{M) and get, using the 
symmetry of Ag 



Notice that Jj^ Nu^hd^g < ||A/'||i:,oc ||m£||2;,2 ||/i||j;,2, so that the right hand side of 
(15) tends to zero. Since h was arbitrary and H'^{M,fi.g) is dense in L'^{M,iJ,g) 
{M being compact), we end up with (understood in weak sense) 



But since ^ > and the spectrum of —Ag + {£, + V) is nonnegative, the unique 
solution is w — Mo = 0, hence w = Uq. The fact that none of this depends the 
selection of / implies convergence in the operator norm. The statement about 
the diffusion distances is now clear from their definition and the fact that the 
kernels are smooth. 




{w — u^)Aghd^g + / + V){'w — u^)hdiJ,g = e / Nu^hdjig. (15) 





-Ag{w - Wo) + (C + V){w - uq) = 0. 



(16) 



